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Abstract 

The  production  of  electricity  in  a  fuel  cell  system  is  associated  with  the  production  of  an  equivalent  amount  of  thennal  energy,  both  for 
large  size  power  plants  and  for  transportation  applications.  The  heat  released  by  the  cells  must  be  removed  by  a  cooling  system,  characterized 
by  its  small  size  and  weight,  which  must  be  able  to  assure  uniform  work  conditions  and  reduce  performance  losses.  Based  upon  realistic 
assumptions,  a  mathematical  model  has  been  developed  to  determine  the  temperature  and  current  density  distribution  in  a  solid  polymer 
electrolyte  fuel  cell  (SPEFC)  stack  as  a  function  of  operating  conditions  and  stack  geometry.  The  model  represents  a  useful  tool  to  identify 
operating  conditions,  such  as  to  have  an  optimal  longitudinal  and  axial  temperature  profile,  so  allowing  the  design  of  cooling  system  and 
bipolar  plates.  In  this  paper,  the  model  has  been  applied  to  determine  the  temperature  profile  of  an  experimental  SPEFC  stack.  The  model  is 
validated  by  comparing  model  results  with  experimental  measurements;  simulated  and  experimental  results  agree  satisfactorily. 
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1.  Introduction 

A  fuel  cell  stack  usually  consists  of  a  number  of  single 
cells  connected  in  electrical  series.  Each  cell,  whose  basic 
elements  are  a  matrix  containing  electrolyte  and  porous  elec* 
trodes  (anode  and  cathode),  is  connected  to  the  adjacent  one 
through  bipolar  plates  ribbed  with  channels  that  allow  the 
flow  of  gases.  Thus,  the  bipolar  plates  are  gas  and  current 
collectors,  and  they  provide  the  separation  between  the  fuel 
and  oxidant  flows  in  adjacent  cells. 

Since  the  electrical  energy  produced  by  the  cells  involves 
a  corresponding  amount  of  thermal  energy,  it  is  necessary  to 
study  an  appropriate  cooling  system  that  can  remove  this  heat 
and  maintain  the  working  temperature  at  a  constant  value.  In 
fact,  high  temperature  gradients  affect  the  cell  performance 
[1].  In  particular,  in  solid  polymer  electrolyte  fuel  cells 
( SFRFCs) ,  temperature  gradients  cause  problems  with  water 
management  and  electrochemical  reactions,  then  reducing  the 
cell  lifetime  [2].  Heat  removal  is  generally  accomplished 
through  the  use  of  cooling  plates,  which  are  usually  placed 
in  every  repeating  unit  of  a  given  number  of  working  plates. 
It  is  evident  that  different  cooling  fluids  (air,  oil,  water,  etc.) 
and  geometric  configurations  (channels,  tubes,  solid  sur¬ 
faces,  etc.)  with  sizes,  shapes  and  materials  as  variables  can 
be  proposed  for  an  efficient  cooling  system.  Often  this  choice 
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is  dictated  by  practical  and  economic  requirements;  a  sim¬ 
plified  cooling  system  is  preferable  in  comparison  with 
solutions  which  farther  increase  the  system  complexity  and 
management. 

Even  if  experimental  analyses  and/or  theoretical  treat¬ 
ments  of  these  problems  have  been  carried  out  by  other 
research  groups,  see  Refs.  [3~6],  only  few  literature  works 
[7,8]  are  devoted  to  the  application  of  these  studies  to  the 
case  of  SPEFC.  In  fact,  this  type  of  fuel  cell  has  been  exten¬ 
sively  studied  mainly  from  the  experimental  viewpoint,  while 
the  theoretical  studies  concerned  exclusively  micromodels 
for  the  single  cell  and  its  components  [9-1 1  ] . 

This  paper  presents  a  three-dimensional  mathematical 
model  applied  to  study  the  current  and  temperature  distribu¬ 
tion  of  an  SPEFC  stack,  whose  cooling  system  is  represented 
by  U-tubes  supported  by  each  working  plate;  the  cooling  fluid 
is  water.  The  model  has  been  applied  to  verify  the  results  of 
an  experimental  stack  built  by  De  Nora  Permelec  SpA  (DNP) 
Milan,  Italy  [  12] .  For  this  purpose,  the  simulation  model  has 
been  studied  appropriately  to  comply  with  the  technical  char¬ 
acteristics  of  the  DNP  stack.  The  results  obtained  by  the 
model  are  in  good  agreement  with  experimental  measure¬ 
ments  and  demonstrate  the  validation  of  the  model  calcula¬ 
tions.  A  sensitivity  analysis  of  the  calculated  temperature 
profiles  as  a  function  of  the  operating  conditions  ami  the 
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cooling  system  geometry  allowed  the  optimization  of  the 
stack  design.  In  particular,  the  cooling  system  warrants  an 
average  temperature  close  to  the  typical  values  for  this  kind 
of  fuel  cell  (70-75  °C)  and  a  uniform  temperature  profile  in 
the  stack. 


2.  Experimental 

An  SPEFC  stack  equipped  with  temperature  and  pressure 
probes  was  set  up  to  allow  monitoring  of  the  stack  operating 
parameters.  Temperature  measuring  devices  (F  thermocou¬ 
ples)  were  also  inserted  in  the  core  of  each  bipolar  plate  to 
control  the  temperature  al  the  centre  of  the  active  cell  area. 
An  aluminium  cooling  plate,  coupled  with  ribbed  graphite 
current  collectors,  constituted  the  bipolar  plate.  Fifteen  cells 
were  assembled  in  one  unit  in  a  standard  filter-press  config¬ 
uration  by  means  of  tie  bolts,  in  order  to  avoid  gas  leakage 
and  to  ensure  proper  electrical  contact.  Humidification  of  the 
reactant  gases  was  controlled  by  the  use  of  external  humidi¬ 
fiers  allowing  gas  saturation  at  the  prefixed  temperature  and 
flow  rates.  Nation®  1 17  membranes  and  DNP  electrodes,  with 
1  mg/cm2  Pt  loading,  were  used  as  the  active  components  in 
each  cell. 

A  schematic  of  the  SPEFC  stack  is  given  in  Fig.  1,  where 
the  repeating  section  of  the  stack  —  cooling  plate,  anodic 
current  collector,  anode,  membrane,  cathode  and  cathodic 
current  collector — constitutes  the  elementary  cell. 


a)  123  cooling  tubes 

-r-,  ..  ...  •)  fS  membrane 

w)  6sl  cooling  plate  collector 
v  I)  □  cathode 

c)  L_J  anodic  cunont  colloclor 

g)|  |  cathodic  current  collector 

d)  H  anode 

Fig.  1.  Schematic  view  of  the  elementary  unit  of  the  SPEFC  stack  for  the 
temperature  distribution  three-dimensional  model. 


3.  Modelling  development 


To  describe  the  distribution  of  current  density  and  temper¬ 
ature  in  an  SPEFC  stack  through  a  simulation  model  some 
assumptions  have  been  made  to  simplify  the  mathematical 
treatment.  All  the  assumptions,  even  if  they  do  not  claim  to 
have  a  general  validity,  have  been  conceived  so  that  they  are 
coherent  with  the  case  considered  and  result  in  a  realistic  and 
appropriate  description  of  the  examined  problem.  The  main 
assumptions  are  summarized  as  follows:  (i)  to  simulate  the 
U-tubes,  the  cooling  system  has  been  considered  as  a 
sequence  of  separate  tubes  of  a  cooling  plate  with  water  flow 
in  alternating  direction  and  perpendicular  to  the  air  flow;  (ii) 
the  pressure  drop  (using  relationships  reported  in  literature 
[  13,14],  a  value  of  841  Pa  has  been  evaluated)  andtheeffects 
associated  with  the  curvature  of  the  U-tubes  have  been 
neglected  ( except  for  some  boundary  conditions  on  the  water 
temperature  introduced  to  simulate  the  continuity  between 
adjacent  tubes);  (iii)  the  edges  of  the  stack  operate  adiabat- 
ically;  (iv)  the  effective  thermal  conductivities  of  the  plates 
in  the  air  flow  direction,  and  in  the  direction  perpendicular  to 
the  air  flow,  are  considered  to  be  equal,  while  the  thermal 
conductivity  of  the  cooling  tubes  is  considered  negligible; 
( v)  stack  symmetry  allows  us  to  consider  only  the  working 
plates  between  two  different  cooling  plates  to  evaluate  current 
and  temperature  distribution;  and  (vi)  temperature  gradients 
in  the  flow  direction  of  the  fuel  are  assumed  to  be  negligible, 
i.e.  the  fuel  heat  capacity  is  much  less  than  the  oxidant  heat 
capacity. 

Based  on  these  assumptions,  the  general  equations  describ¬ 
ing  the  energetic  balances  relative  to  the  cooling  plate,  the 
working  plates,  the  cooling  water  and  the  process  air  are  given 
by: 

(i)  cooling  plate 
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(iv)  process  air 
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with  the  boundary  conditions 
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ar  _Q  ar  _Q 
l,-o  ’  ^  1,-1, 

(adiabatic  conditions) 


t*=ta+ 2tc* 

and 


(5a,  b) 
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(5c,  d) 

(adiabatic  conditions) 

(6) 

(symmetric  condition) 

r«l,.o.,.o=rWJn 

(7) 

Twlx  =  fc,..V=/£y/WC  =  Tw|Jr  =  o,.V  =  t/-  Dty/lVC 

(8a) 

1  £7  £  NC  j  odd 

1  jr-o.  y  -JL^NC  =  Fwl,-o,.v- (/- l)£,/WC 

(8b) 

1  £  /"  £  NC  j  even 

Tpljr-0-rpj. 

(9) 

where 

hjj  ~ 


£0.001 


Then,  the  Eqs.  (l)-(4),  with  the  boundary  conditions 
(5)-(9),  are  solved  to  determine  the  temperature  distribu¬ 
tions  for  the  cooling  plate,  working  plates,  cooling  water  and 
process  air.  The  procedure  is  iterative:  the  final  temperature 
distribution  is  obtained  when  the  differences  between  the 
values  of  temperature  calculated  for  an  iteration  loop  and 
those  of  the  previous  loop  (an  initial  approximation  is  used 
for  the  first  loop)  are  lower  than  a  specified  tolerance,  in 
accordance  with  the  following  criterion  of  convergence: 


I  'f'r.ij  ^t,ij 

\TklJ+TX)kJJ 


£0.002 


The  finite-difference  method  replaces  the  differential  Eqs. 
(l)-(4)  with  algebraic  equations,  substituting  the  deriva¬ 
tives  with  local  differences  evaluated  in  the  grid  points.  By 
using  the  indices  i,j,  k  to  denote  the  subdivisions  along*-,  y- 
and  z-axes,  which  respectively  represent  the  direction  per¬ 
pendicular  to  the  air  flow,  the  process  air  direction,  and  the 
stack  (axial)  direction,  the  finite-difference  equations  are: 

(i)  cooling  plate  (fc=  1 ) 


<“  =  (.+ »c + tm  4-  +  iccB/2 

are  the  total  thicknesses  of  the  cooling  and  working  plate, 
respectively 5  k  is  the  effective  thermal  conductivity,  C  is  the 
heat  capacity,  m  is  the  mass  flow  and  h  is  the  heat-transfer 
coefficient  (for  the  other  definitions  the  reader  can  refer  to 
Section  6) 

The  ‘plus-or-minus’  terms  appearing  in  Eqs.  (1)  and  (3) 
take  into  account  the  alternating  direction  of  the  water  in  the 
U-tubes;  in  fact,  in  accordance  with  assumption  (i)  in  this 
Section,  the  water  flows  perpendicularly  to  the  air,  going  from 
left-side  to  right-side  of  the  plates  in  odd  tubes  (sign  - )  and 
in  the  opposite  direction  in  even  tubes  (sign  + ) .  The  bound¬ 
ary  conditions  (8a)  and  (8b)  are  also  valid  for  odd  and  even 
tubes,  respectively  ( NC  is  the  number  of  cooling  tubes) . 

The  procedure  used  to  solve  these  equation1:  is  based  on  a 
finite-difference  method  that  subdivides  each  plate  in  a  num¬ 
ber  of  grids,  which  are  each  small  enough  so  that  temperature 
and  concentrations  variations  may  be  neglected.  Once  the 
average  current  density  is  specified,  the  model  calculates  the 
voltage  of  each  working  plate.  This  calculation  utilizes  gen¬ 
eral  literature  relationships  ( 15],  opportunely  adapted  to  the 
SPEFC,  through  empirical  correlations  that  give  the  different 
overpotential  losses  (activation,  diffusional,  ohmic)  as  a 
function  of  the  electrode  characteristics. 

The  current  density  distribution  on  the  working  plates  is 
determined  as  a  function  of  the  local  (on  the  grids)  reactant 
compositions.  A  solution  that  satisfies  the  input  specifications 
(average  current  density,  reactant  utilizations  and  flow  rates, 
etc.)  is  obtained  when  the  differences  belveen  the  local  val¬ 
ues  and  an  initial  approximation  (I0kjj)  of  the  current 
density  are  within  a  given  tolerance: 


A\TklJ_ |  A,r*^+|  C'Tk+ 

~~  CfTwf  +  i.ij + D\  ( TwtJ — TWilJ) 

-El(rkJ_u+r^+lJ)= o  do) 

( ii)  working  plates  (2  £  k^NK  + 1 ) 

— A2  ,  +  B2T^,j—A2TkjJ+ ,  —B3Tk-Uj 
~  BtTk+  tiJ + D2(  Tp*_  —  T^_  ,) 

(11) 

(iii)  cooling  water 

TwU=TWUJ+  (7^-7^)  [  1  -cxp(  -  (12) 

(iv)  process  air  (2  £fc£AflC+l) 

^P*jj  —  1  "**  ( +  U J  ~  ^Vmj-  » ) 

X[l~exp(-0PW)]  (13) 

where  TkjJ,  TWMand  T9kiy  with  1  £Jt£AttT+ 1, 1  £i£  JVXand 
1  £  j  £  NY  represent  the  temperatures  of  the  {Mates  (k  —  1  for 
cooling  plate),  of  the  cooling  water  ami  of  the  process  air  in 
the  grid  points,  NK  being  the  number  of  cells  between  two 
different  cooling  plates  and  NX,  NY  the  number  of  subdivi¬ 
sions  along  *  and  y,  respectively,  and 

T  =/  T-‘-J  fOT  J  0dd 
l Tm.u  for;  even 

The  above  coefficients  are  defined  as  follows: 
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Tabic  1 

Coefficients  for  the  finite-difference  Eq.  ( 1 1 ) 


Coefficients  *~2,  if  NK-  l  k=2.i(NK>l  2  <k£NK  k™NK+\ 
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B*  0  C2  C2  0 


When  k*NK  +  1  the  term  fi47'*+ \jj  is  replaced  with  Ci  TUr 


A _*££ 

1  Ay2 

(14a) 

j"a/ 

(14b) 

Bt  =  C,‘+Cf  +  E, 

(15) 

X-A.U 

a=— — 

1  AsA'u 

(16a) 

1.  AX 

Q  ArAL 

(I6b) 

2  Az"x  A<-B-u 

(16c) 

C£|. 

'  PJLx 

(17a) 

Dj=S & 

PpDy 

(17b) 

„ 

E'  Ajt 

(18a) 

5ft 

II 

(18b) 

E3  =  2A,+2£i 

(18c) 

E,  =  2A2  +  2E2 

(I8d) 

Ouj-IVT-Vujilui 

(19) 

y'v,J  fM  r 

WwW 

(20a) 

™pcp 

(20b) 

The  coefficients  B2t  £3  and  fi4,  which  appear  in  Eq.  ( 1 1 ) 
for  a  generic  working  plate  {2£k£NK+ 1),  are  listed  in 
Table  1.  They  depend  on  the  thermal  conductivities  of  adja¬ 
cent  plates,  and  therefore  vary  with  the  location  of  the  con¬ 
sidered  plate.  In  fact,  for  the  first  and  last  working  plates 
(k  —  2  and  k**NK+ 1 ),  the  conductivity  between  the  cooling 
plate  and  the  anodic  ( CV)  or  cathodic  section  of  the  working 
plate  ( Cic)  must  be  considered.  For  the  other  working  plates 
(2&kgNK),  however,  only  the  conductivities  between 
working  plates  (C2)  come  into  play. 

The  boundary  conditions  (5)-(9)  become 


(21a) 

TkJJ*Y+ 1  “  Tkj jvy 

(21b) 

Tk.oj=TUJ 

(21c) 

(21d) 

Tnk+2JJ  =  Tnk+  1  ,ij 

(22) 

rwol  -  rwin 

(23) 

T„0J=TW)J^  (for;  even) 

(24a) 

U  =  W  »  <  f°r  j  “W) 

(24b) 

^,0-y-p.in 

(25) 

By  considering  the  number  of  subdivisions  in  the  air 
flow  direction  equal  to  the  number  of  cooling  tubes 
( 1  £j^NY=NC),  it  has  been  possible  to  express  the  conti¬ 
nuity  conditions  between  adjacent  tubes  (24a)  and  (24b), 
which  require  that  the  inlet  water  temperature  in  each  tube 
(or  division)  is  equal  to  the  exit  water  temperature  of  the 
previous  tube  (division). 

The  simultaneous  solution  of  Eqs.  ( 10)— ( 13)  is  complex, 
because  the  number  of  unknown  temperatures  flu*  arK* 
T^j  is  equal  to  NX- NY-  (NK  + 1 ) ,  which  is  also  the  number 
of  equations  to  be  solved.  To  reduce  the  storage  of  matrix 
coefficients  and  the  computational  time,  the  computer  code 
has  been  arranged  to  solve  the  Eqs.  ( 10)— ( 13)  for  a  fixed 
subdivision  i  (that  corresponds  to  a  slice  in  the  direction  of 
air  flow).  In  practice,  the  three-dimensional  problem  can  be 
reconsidered  as  a  two-dimensional  one,  solving  only 
NX-  {NK  + 1 )  equations  for  each  division  NY. ,  rather  than 
NX •  AT-  ( NK  + 1 )  simultaneous  equations.  For  this  purpose, 
at  the  first  calculation  loop  the  thermal  conductivities  of  the 
plates  in  the  direction  perpendicular  to  the  air  flow  are 
neglected  (kxA  =*  k*  =  0,  that  is  E,  =  £2  =  0  in  Eqs.  ( 10)  and 
(11))  and  a  first  approximation  of  the  water  temperatures 
can  be  calculated  from  specified  (experimental)  inlet  and 
outlet  water  temperatures  in  the  U-tube,  assuming  a  linear 
variation  with  the  length  of  the  tube.  Beginning  with  the 
second  loop,  the  conductivities  are  taken  into  account 
(kxA = kyA  and  kxB = kyB  according  to  assumption  (iv)  in  this 
Section)  ;  besides,  for  each  division  i,  the  unknown  temper¬ 
ature  of  the  t  —  1  divisions  is  replaced  with  the  calculated 
value  (from  the  previous  division)  and  the  i 4- 1  temperature 
is  replaced  with  the  corresponding  value  calculated  from  the 
previous  loop. 


G.  Maggio  ex  al  /  Journal  of  Power  S  crees  62  (1996)  167-174 


171 


TTic  effects  associated  with  the  thermal  conductivities  of 
both  working  and  cooling  plates,  the  conduction  between 
adjacent  plates,  and  the  separate  cooling  effects  of  the  process 
air  and  cooling  water  have  been  examined.  In  particular,  the 
effective  thermal  conductivities  of  the  plates  have  been  cal¬ 
culated  as  a  function  of  both  the  thickness  and  the  conductiv¬ 
ity  of  the  individual  elements  of  the  plate  (anodic  and 
cathodic  current  collectors,  anode,  cathode  and  membrane). 
Tbe  possibility  to  consider  dissymmetries  between  anodic 
and  cathodic  elements,  in  terms  of  thickness  and  thermal 
conductivity  (different  materials),  has  also  been  introduced 
in  the  model. 

As  regards  to  the  direction  of  the  reactant  gas  flows  the 
model  allows  various  choices,  corresponding  to  differentcon- 
figurations  (perpendicular,  parallel  equi-current  or  counter- 
current),  that  do  not  involve  modifications  of  the  assumption 
(vi)  of  this  Section. 


4.  Results  and  discussion 

The  model  has  been  applied  to  verify  the  distribution  of 
temperature  in  an  experimental  SPEFC  stack.  The  stack  con¬ 
sists  of  IS  cells  with  a  surface  area  of  225  cm2,  working  at  a 
current  of  100  A.  Each  cell  is  cooled  (in  this  case  the  number 
of  working  plates  between  two  cooling  plates  is  NK—  1)  by 
a  U-tube  with  10  passages  (simulated  as  single  tubes),  that 
have  a  cylindrical  geometry  and  an  internal  diameter  of  6 
mm;  the  total  water  flow  is  0.3  Nm3/h  (20  Nl/h  per  cell) 
and  the  inlet  water  temperature  is  7 1 .3  °C.  The  other  operating 
conditions  and  geometric  features  of  the  examined  stack,  used 
by  the  model  as  input  data,  are  listed  in  Table  2.  The  elec¬ 
trodes  are  considered  to  have  the  same  characteristics  as 
CNR-TAE  electrodes  [  16] . 

The  effective  thermal  conductivity  between  cooling  and 
working  plate  ( anodic  or  cathodic  section) ,  calculated  by  the 
thicknesses  and  the  conductivities  of  the  individual  elements, 
is  1 5.7  W/  ( m  K) .  This  is  due  to  the  fact  that  no  dissymmetries 
between  anodic  and  cathodic  elements  have  been  considered 
(same  thickness)  and  the  material  of  the  cooling  plates  is  the 
same  the  working  plates.  Besides,  the  working  and  the 
cooling  plates  effective  thermal  conductivities  (in  air  or  per¬ 
pendicular  direction)  were  136.5  and  205.7  W/(m  K), 
respectively;  in  particular,  the  last  one  is  equal  to  the  value 
relative  to  the  cooling  plate  collector  because  the  thermal 
conductivity  of  the  tubes  was  neglected.  The  effective  thermal 
conductivity  between  two  working  plates,  in  this  case,  has 
not  been  considered,  since  the  working  plate  is  in  contact  with 
a  cooling  plate  both  in  the  upper  and  in  the  lower  side  (Fig. 
1)  and  there  is  no  current  collector  between  working  plates 

(fccB~0). 

The  value  of  the  reversible  potential  and  the  overpotentials 
calculated  by  the  model  as  a  function  of  the  electrode  char¬ 
acteristics  are  reported  in  Table  3.  The  limiting  current  asso¬ 
ciated  with  the  concentration  overpotential  resulted  equal  to 
1.17  A/cm2.  The  anodic  activation  overpotential  is,  gener- 


Table  2 

Model  input  data 


Parameter  description 

Value 

Working  plate  size 

15  cmx  15  cm 

Prefixed  average  current  density 

444  mA/cm2 

Oxygen  utilization  coefficient 

69.6% 

Hydrogen  utilization  coefficient 

66.3% 

Cathodic  total  pressure 

3.0  atm 

Anodic  total  pressure 

2.5  am* 

Number  of  cooling  tubes 

10 

Fuel  inlet  temperature 

68  °C 

Cathode  humidification  temperature 

71  °C 

Anode  humidification  temperature 

70  °C 

Inlet  cooling  water  temperature 

71.3  °C 

Inlet  process  air  temperature 

80  CC 

Total  cell  number 

15 

Finite-differences  along  y 

12 

Oxidant 

Air 

Cooling  tubes  external  diameter 

8  mm 

Electrode  thickness 

0.3  ram 

Membrane  thickness 

0  !7  mm 

Cooling  plate  collector  thickness 

2mm 

Current  collector  thickness 

1.5  mm 

Catalyit  weight 

1  mg/cm2 

Catalyst  utilization 

30% 

Catalyst  surface  area 

600  cnr7mg 

Total  water  flow 

0.3  Nra3/h 

Electrode  therm’d  conductivity 

15.56  W/(m  K) 

Membra**  thermal  conductivity 

0.21  W/(mK) 

Coojirii  plate  collector  thermal  conductivity 

205.7  W/(m  K) 

'Current  collector  thermal  conductivity 

30W/(mK) 

Working  plate  density 

2.65  g/cm3 

Cooling  plate  density 

2.65  g/cm3 

Working  plate  heat  capacity 

0.405  kcal/ (kg  *’C) 

Cooling  plate  heal  capacity 

0.405  kcal/fkg’C) 

Table  3 

Electrode  overpotentials  calculated  by  the  model 

Parameter 

Value  (mV) 

Reversible  potential 

1169.8 

Cathodic  activation  overpotential 

472.2 

Anodic  activation  overpotential 

0.0 

Ohmic  overpotential 

103.6 

Concentration  overpoteniial 

3.6 

ally,  a  function  of  the  carbon  monoxide  concentration  and 
can  be  obtained  by  empirical  equations,  hi  the  SPEFC  case, 
the  fuel  stream,  going  into  the  stack,  must  not  contain  carbon 
monoxide  (CO  acts  as  a  catalyst  poison),  then  we  assumed 
the  anodic  activation  overpotential  to  be  zero.  A  first  approx¬ 
imation  of  the  cell  voltage  (590.4  mV)  is  calculated  by  the 
model  using  the  values  of  the  parameters  listed  in  Table  3. 
Then,  this  value  is  revised  to  obtain  the  convergence  of  tbe 
calculated  current  densities  at  the  prefixed  average  current 
density  (444  mA/cm2),  resulting  in  a  final  cell  voltage  equal 
to  604.6  mV. 

The  model  validation  has  been  accomplished  comparing 
the  cal  .ulated  average  cell  temperature,  cell  voltage  and  out¬ 
let  water  temperature  with  the  measured  values  (Table  4). 
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Table  4 

Comparison  between  model  and  experimental  results 


Result 

Experimental 

Calculated 

Difference 

Average  temperature  (°C) 

75.3 

74.69 

-0.81% 

Average  cell  voltage  (mV) 

587.8  * 

604.6 

+  2.86% 

605.  i  6 

-0.08% 

Outlet  water  temperature  (°C) 

73.9 

74.3 

+0.54% 

*  Average  voltage  calculated  on  15  cells. 
b  Average  voltage  calculated  on  the  first  14  cells. 

From  Table  4  it  is  evident  that  the  agreement  between  the 
simulation  and  the  experimental  results  is  very  satisfactory. 
The  differences  between  experimental  and  calculated  values 
of  the  cell  temperature  and  the  outlet  water  temperature  are 
very  low  ( less  than  1%).  The  difference  between  the  average 
cell  voltage  calculated  by  the  model  and  the  experimental 
one  (the  experimental  voltages  of  every  individual  cell  are 
reported  in  Table  5)  is  a  little  higher  (about  3%),  but  this 
disagreement  is  most  likely  due  to  the  poor  performances  of 
the  last  cell  (whose  voltage  is  346  mV,  due  to  an  inadequate 
water  balance  of  the  cell),  that  for  its  accidental  occurrence 


Table  5 

Experimental  cell  voltages 


Cell  number 

Cell  voltage  (mV) 

1 

630 

2 

588 

3 

610 

4 

616 

5 

612 

6 

590 

7 

618 

8 

613 

9 

591 

10 

603 

II 

605 

12 

611 

13 

603 

14 

581 

15 

346 

Coaling  water  temperature  (C) 


H20  out 


Fig.  2.  Cooling  water  temperature  distribution  forthe  SPEFC  stack  calculated 
by  the  model  (input  data  listed  in  Table  2). 


was  not  foreseeable  by  the  model.  In  fact,  when  the  average 
experimental  voltage  is  calculated  on  the  other  cells,  the  dif¬ 
ference  reduces  to  0.5  mV  (less  than  0.1%). 

The  calculated  temperature  and  current  density  distribu¬ 
tions  are  represented  in  Figs.  2-5.  In  particular,  the  temper¬ 
ature  profile  of  the  cooling  water  and  that  of  the  cooling  plate 
are  shown  in  Figs.  2  and  3,  respectively;  the  temperature 
profile  of  the  working  plate  and  its  current  density  are 
depicted  in  Figs.  4  and  5,  respectively.  In  accordance  with 
the  design  of  the  DNP  experimental  stack,  the  flows  of  the 
process  gases  have  been  considered  to  be  parallel equi-current 
and  their  inlet  side  is  at  the  back  of  the  plates.  The  number  of 
grids  considered  for  the  finite-difference  method  is  equal  to 


Cooling  plate  temperature  (C) 


H20  out 


Fig.  3.  Cooling  plate  temperature  distribution  for  the  SPEFC  stack  calculated 
by  the  model  (input  data  listed  in  Table  2). 


Working  plate  temperature  (€) 


H20  out 


Fig.  4.  Working  plate  temperature  distribution  for  the  SPEFC  stack  calcu¬ 
lated  by  the  model  (input  data  listed  in  Table  2). 
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Current  density  (mA/cm2) 


H20  out 


Fig.  5.  Current  density  distribution  for  the  SPEFC  stack  calculated  by  the 
model  (input  data  listed  in  Table  2). 

12X10  (the  second  number  represents  the  subdivisions 
along  air  flow  direction  and  is  equal  to  the  number  of  cooling 
tubes),  because  this  subdivision  proved  sufficient  for  an  accu¬ 
rate  solution. 

The  average  temperatures  of  the  cooling  and  working 
plates  are  73.94  and  74.69  °C,  respectively.  Besides,  the  cor¬ 
responding  average  temperature  gradients  (about  0.70  and 
0.53  °C ) ,  demonstrate  that  a  uniform  temperature  distribution 
is  achieved  in  the  stack.  The  average  temperature  of  the  cool¬ 
ing  water  is  72.81  °C,  while  that  of  the  process  air  is  close  to 
the  average  cell  temperature  (74.69  °C).  The  average  current 
density  is  444  mA/cm2,  and  its  values  range  from  370  to 
505  mA/cm2. 


S.  Conclusions 


In  this  paper  a  three-dimensional  mathematical  model  has 
been  proposed  to  study  the  current  density  and  temperature 
profiles  in  an  SPEFC  stack.  The  model  allows,  through  a 
sensitivity  analysis,  to  study  the  influence  of  some  geomet¬ 
rical  (size,  shape  and  number  of  the  channels  of  the  cooling 
plate,  etc.)  and  working  parameters  (inlet  temperatures  and 
flow  rates  of  the  cooling  water  and  air,  anodic  and  cathodic 
pressure,  etc.)  on  the  stack  performance. 

The  model  validation  has  been  successfully  accomplished 
as  evidenced  by  the  comparison  between  the  calculated 
results  provided  by  the  model  and  the  experimental  ones 
corresponding  to  an  SPEFC  stack  with  an  original  configu¬ 
ration  built  by  De  Nora  Permelec  SpA  in  Milan,  Italy  [  12]. 
The  slack  consists  of  1 5  cells  with  a  surface  of 225  cm2  cooled 
by  the  water  flowing  through  a  U-tube  (total  water  flow  is 
0.3  Nm3/h  at  inlet  temperature  of  71.3  °C)  and  is  designed 
for  a  total  electrical  power  of  1  kW. 

The  average  cell  voltage  experimentally  measured  for  the 
stack  is  about  600  mV  at  100  A,  in  good  agreement  with  the 
value  determined  by  the  model  (604.6  mV).  The  cell  tem¬ 
peratures  calculated  by  the  model  (average  value  is  74.7  °C) 
are  in  agreement  with  the  experimental  values  as  recorded  by 
thermocouples  placed  inside  the  stack,  with  differences  lower 
than  1%  on  the  average  values.  Additionally,  temperature 


profiles  relative  to  the  cooling  plate,  the  process  air  and  the 
cooling  water  have  been  determined  by  the  model  calcula¬ 
tions:  average  values  are  respectively  equal  to 73.9, 74.7  and 
72.8  °C. 

The  choice  of  appropriate  operating  conditions  and  geo¬ 
metrical  features  of  the  stack  allowed  to  obtain  an  average 
temperature  of  the  working  plates  (about  75  °C)  close  to  the 
value  required  for  this  kind  of  fuel  cell.  A  uniform  distribution 
of  the  temperature  in  the  stack  has  also  been  achieved  (the 
gradient  of  temperature  is  less  than  4  °C  both  in  each  plate 
and  in  the  axial  direction  of  the  stack) . 


6.  List  of  symbols 


A1xA2 

Z?1,  fl2t  ^3.  ^4 
C 

CAC,c.C2 

Dj,  Dj 

Eu  E2>Ey,  £4 
F 
h 
/ 
it 

L 

m 

n 

NC 

NK 

NX 

NY 

P 

S 

t 

T 

V* 

V 

x,y,z 


see  Eqs.  ( 14a)  and  (b) 
see  Eq.  ( 15)  and  Table  1 
heat  capacity,  J/  (kg  K) 
see  Eqs.  (16a)-(c) 
see  Eqs.  ( 17a)  and  (b) 
see  Eqs.  (18a)-(c) 

Faraday  constant  (96487  C/g-eq.) 
heat-transfer  coefficient,  J/  ( s  cm2  K) 
cell  current  density,  A /cm2 
effective  thermal  conductivity,  J/ (s  cm 
K) 

slice  length,  cm 
mass  flow,  kg/s 

number  of  electrons  transferred  by  the 
ceil  reaction  (2  g  eq./mol) 
number  of  cooling  tubes 
number  of  cells  between  two  cooling 
plates 

subdivisions  along  jc-axis 
subdivisions  along  y-axis 
pitch  of  tube,  cm 
tube  perimeter,  cm 
plate  thickness,  cm 
temperature,  K 
Aff/nF,  V 
cell  voltage,  V 

perpendicular  to  the  process  air,  process 
air  and  stack  (axial)  direction, 
respectively 


Greek  symbols 

A  H  reaction  enthalpy  (  J/mol ) 

0  see  Eqs.  (20a)  and  (b) 

n  seeEq.  (19) 

Superscripts 

A  cooling  plate 

B  working  plate 

L  lower  plate 

U  upper  plate 
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Subscripts 

a 

c 

cc 

ct 

in 

m 

P 

w 


anode 

cathode 

cuuent  collector 
cooling  tube 
inlet 

membrane 
process  air 
cooling  water 
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